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AN  ALGORITHM  FOR  IMPROVED  GATING 
COMBINATORICS  IN  MULTIPLE-TARGET  TRACKING 


I.  Introduction 

In  this  paper  we  describe  a  method  for  significantly  reducing  the  computational  com¬ 
plexity  required  for  observation-track  gating  in  multiple  target  tracking.  We  define 
the  gating  process  as  follows:  given  a  set  of  observations  and  Ht  jtracks,  identify 
all  observation-track  pairs  whose  scores  fall  above  a  chosen  threshold.  The  score  for 
observation-track  pair  (t,j)  is  defined  as  the  function; 


s.,(dx.„r,) 


exp(-dX/,r-^dX.y/2) 

(2T)^yir;j 


(1) 


where  d  is  the  measurement  dimension,  Fj  is  the  residual  covariance  matrix  of  the  track 
j,  and  dXij  is  the  residual  vector  of  the  pair  (t.ji),  and  these  arguments  are  taken  to 
be  valid  at  the  same  time.  In  addition,  throughout  this  study  we  shall  use  the  score 
threshold,  Sj,  chosen  for  every  observation  j  such  that: 


Sii  <  S,  . 


(2) 


In  part  because  Eq.  (1)  is  very  expensive  to  evaluate,  much  of  the  previous  work  on 
gating  has  emphasized  the  use  of  intermediate  or  “coarse”  gating  criteria  which  replaces 
the  calculation  of  Eq.  (1)  with  a  function  that  is  computationally  cheaper  to  evaluate. 
The  result  is  the  identification  of  a  superset  of  candidate  pairs  which  includes  the  pairs 
that  satisfy  Eq.  (2).  We  denote  the  subset  satisfying  Eq.  (2)  as  either  the  pairs  which 
correlate  at  gating  or  the  pairs  which  satisfy  the  final  gate.  Typically  this  pre-processing 
includes  defining  a  gating  volume,  Vq,  based  on  the  numerator  of  the  right-hand  side 
(RHS)  of  Eq.  (1)  so  that  a  coarse  correlation  measure  is  evaluated  assuming  a  Gaussian 
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distribution  as  in  Eq.  (1).  For  example,  the  pairs  pass  one  gate  if  "yVc  < 
where  7VG  is  a  threshold  that  may  be  obtained  from  a  table  or  an  error  function  inte¬ 
gration  for  a  certain  probability  of  correlation  [1,  page  97]  [2,  pages  88ff].  These  pairs 
might  also  be  pre-processed  by  coarser  gating  criteria  with  larger  Vc  but  which  are 
cheaper  to  evaluate  [2,  page  91],  [1,  page  97],  Ideally,  one  completes  the  gating  calcula¬ 
tion  by  computing  5,j  for  the  set  of  candidate  pairs  and  performing  the  comparison  in 
Eq.  (2).  Overall  processing  work  would  be  reduced  since  only  the  pairs  with  sufficiently 
high  coarse  correlation  values  would  be  re-evaluated  using  the  numerically  expensive 
function  of  Eq.  (1). 

These  techniques  address  the  coefficient  of  the  scaling  but  not  the  scaling  itself  since  they 
explicitly  apply  a  coarse  correlation  function  to  all  NtNh  possible  pairs.  If  the  number 
of  pairs  explicitly  evaluated  by  coarse  gating  are  of  order  (NtNr),  then  for  sufficiently 
large  Nt  and  Nr,  real-time  processing  can  be  precluded  on  any  computer  even  with 
the  use  of  numerically  simple  coarse  correlation  functions.  More  recently  Nagarajan  et 
al.[l,  pages  97-98]  have  recognized  that  one  can  can  reduce  the  number  of  times  that  a 
coarse  correlation  function  is  evaluated.  An  algorithm  is  presented  for  doing  so,  but  the 
analysio  of  scaling  is  incomplete.  In  particular,  neither  cost  nor  scaling  equations  of  the 
gating  process  are  presented  in  terms  of  parameters  important  in  multitarget  tracking, 
e.g.  Nt  or  Nr.  Partitional  spatial  clustering  methods  have  been  proposed  to  efficiently 
perform  gating  by  exploiting  features  of  the  target  distribution[3]'.  However,  improved 
scaling  cannot  be  guaranteed  for  these  data  dependent  algorithms.  The  objective  of  this 
paper  is  to  present  an  efficient  gating  algorithm  and  to  analyse  its  scaling.  We  shall  show 
that  the  overall  algorithm  scales  significantly  better  than  quadratic  even  when  reports 
have  unequal  timestamps  within  a  scan.  Our  algorithm  is  compatible  with  virtually 
all  of  the  previous  work  on  auxiliary  gating  criteria  and  coarse  gates.  As  a  secondary 
objective  we  present  an  inexpensive  algorithm  for  calculating  a  coarse  time  independent 
gate  volume  for  the  criterion  of  Eq  (2). 


II. Preliminary  Details 

By  a  d-dimensional  report,  or  observation,  we  mean  a  set  of  d  elements  measured  simul¬ 
taneously  at  some  specified  time.  We  call  this  time  the  validity  time  of  the  observation 
or  the  timestamp  of  the  observation.  We  require  that  the  timestamps  fall  within  a  period 
of  time  called  a  scan  of  length  r,  where  the  times  of  the  reception  of  the  first  and  the 
last  observations  fix  the  beginning  and  the  end  of  the  scan.  A  track  is  an  estimate  that 
in  some  sense  converges  to  the  “true  trajectory”  as  the  number  of  observations  correctly 
matched  to  the  track  increases.  There  are  I  position  components  to  each  track  at  any 
one  time.  Similarly,  the  observations  are  ^-dimensional  in  position,  with  I  <  d. 

In  our  analysis  and  numerical  simulations  we  employ  algorithms  which  find  near  neigh¬ 
bors  of  points  in  ^-dimensional  position  space,  where  nearness  is  defined  by  the  Eu¬ 
clidean  metric.  They  are  used,  for  example,  to  find  the  tracks  with  mean  positions  near 
a  given  observation  position.  Essential  is  that:  (1)  these  algorithms  find  all  the  neigh¬ 
bors  in  some  expected  optimal  or  near  optimal  time,  and  (2),  that  their  performance  be 
relatively  insensitive  to  spatial  distributions.  After  examining  several  search  algorithms, 

'For  cluster  type  distinctions,  see  [4] 
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we  selected  a  BLD-enhanced  fc-d  search  tree  for  our  test^J^][6].  These  data  structures 
are  known  to  have  a  worst  case  single-query  search  time  scalmg  of  0(N^~‘  -|-  Nc)i  where 
N  is  the  size  of  the  database  being  searched  and  No  is  the  number  of  near  neighbors 
returned.  In  many  applications  (including  our  tests)  their  average  case  search  time 
scaling  approaches  O(log  +  No)-  For  a  background  on  multidimensional  search  struc¬ 
tures  the  reader  might  consult  Bentley  [7],  Friedman  et  al.  [8]  and  [9],  Hanan  [10]  and 
Mehlhorn[ll].  For  an  introduction  to  these  algorithms  the  reader  might  consult  [12]. 


III.  Problem  Overview  in  Light  of  Combinatorics 

Given  a  set  of  Nj  tracks  and  a  set  of  Nn  reports,  there  are  at  most  NxNh  scores 
Si]  which  can  be  formed.  Of  these,  a  fraction  q  of  them  will  fall  above  the  thresholds 
and  satisfy  Eq.  (2),  where  q  could  be  as  low  as  l/Nx  or  I/Nr  or  smaller.  Ideally  we 
would  only  calculate  the  qNxNji  scores;  at  worst  we  would  calculate  all  NjNfi  scores. 
An  example  of  a  quadratic  scaling  approach,  for  the  case  when  reports  have  unequal 
timestamps,  is  the  following  technique:  integrate  the  equations  of  motion  of  each  of  the 
Nx  tracks  to  the  times  of  each  of  the  Nr  reports  and  compute  the  scores.  For  each 
report,  keep  those  scores  that  are  above  the  desired  threshold.  The  dominant  cost  of  this 
is  the  0{NrNx)  score  calculations  and  integrations.  Of  course,  if  each  score  calculation 
is  replaced  by  a  coarse  gate  calculation,  the  scaling  is  still  quadratic. 

Intuitively,  when  tracks  and  reports  ai  2  at  the  same  time,  tracks  and  reports  that  are 
close  in  position  tend  to  be  correlated.  Of  course,  Eq.  (1)  specifies  the  meaning  of 
correlation  at  gating  and  shows  that  other  parameters  in  addition  to  mean  positions 
must  be  considered.  The  covariances  will  in  part  determine  a  gate  volume  around  the 
mean  positions,  and  we  conceptually  relate  gating  to  geometry  by  saying  that  reports 
and  tracks  which  gate  with  each  other  are  those  pairs  with  intersecting  gate  volumes. 
Let  Nn  be  the  number  of  tracks  per  report  that  should  gate,  as  determined  by  Eqs.  (1) 
and  (2).  Let  the  gating  volumes  be  determined  ideally  in  the  sense  that  the  set  of  pairs 
which  should  gate  by  Eqs.  (1)  and  (2)  is  identical  to  the  set  of  pairs  which  gate.  Let  p 
be  the  object  density  and  Vjo  the  ideal  gating  volume  per  report.  Let  the  average  of  a 
quantity  X  over  all  the  reports  be  given  by  X.  Then  the  total  number  of  gating  pairs 
is  N^Nr  =  pV iqNr  =  qNxNR. 

The  prescription  herein  for  calculating  “not  much  more”  than  the  required  qNxNR 
scores  involves  in  part  using  estimates  of  Vjo,  say  Vq,  in  a  search  structure  for  identifying 
the  pairs.  We  have  chosen  a  neighborhood  search  volume  to  be  spherical  in  nature, 
although  it  is  not  necessarily  optimal.  A  search  radius  Ro  specifies  the  search  volume 
Vo-  We  denote  Ro  by  Rq  when  a  given  report  has  the  same  timestamp  as  the  track 
file  to  be  searched.  When  the  number  of  correlations  that  should  be  made  is  small, 
i.e.  when  qNx^R  is  not  comparable  to  NxNr,  then  pV jqNr  is  also  small.  Assume  we 
can  find  an  Rq  per  report  such  that:  (1)  the  actual  search  neighborhood  per  report,VG, 
includes  Via  and  (2)  pV q  is  comparable  to  pV jq. 

When  there  is  a  distribution  in  time  of  tracks  and  reports  throughout  a  scan,  then 
the  required  search  radius  Rq  might  on  average  define  a  search  neighborhood  so  much 
larger  than  Vio  such  that  the  number  of  candidate  pairs  found  is  no  longer  comparable  to 
qNxNR.  After  all,  we  cannot  merely  superimpose  the  tracks  and  reports  and  ask  which 
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error  ellipses  intersect  since  evaluation  of  Eq.  (1)  requires  that  the  function  arguments 
correspond  to  the  same  time.  If  we  instead  ask  which  ellipses  would  intersect  if  they 
were  at  the  same  time,  then  the  gating  volume  around  the  report’s  position  should  take 
into  account  bounds  on  the  location  possibilities  of  the  objects  due  to  their  dynamics 
and  the  time  differences  between  the  report  and  the  tracks.  In  this  case  the  estimate  of 
the  gate  volume  is  zJso  time  dependent,  i.e.  Vq  =  Vc{^T,Ro)-,  and  we  model  its  search 
radius  as 

Rg  =  Rq  +  a|(5T'|  ,  (3) 

where  a  is  some  upper  bound  on  the  velocity  and  ST  is  the  maximum  time  difference 
possible  between  any  track  and  a  report  within  the  scan  r  and  possibly  equal  to  r  itself. 
Thus  scaling  could  depend  on  the  two  parameters  on  the  RHS  pf  Eq.  (3).  In  discussing 
our  algorithm  we  find  that  two  limiting  cases  of  Eq.  (3)  are  particularly  convenient; 


qI^TI  <<  Rq, 

(4a) 

a\ST\  >>  Rq. 

(46) 

We  denote  the  former  case  as  the  “small  scan  length  search”  case  and  the  latter  case 
the  “large  scan  length  search”  case.  Because  of  Eq.  (3)  and  Eq.  (4a),  the  conclusions 
of  the  next  section  on  the  zero  scan  length  case  apply  to  the  small  scan  length  case. 
The  following  section  introduces  the  general  algorithm  through  a  discussion  of  the  large 
sran  length  case.  The  mathematical  analysis  of  the  cost  for  any  value  of  a|^Tl/jRo 
is  presented  afterwards.  We  also  prove  that  the  general  algorithm  has  its  worst  case 
scaling  in  terms  of  a  and  r  when  Eq.  (4b)  holds. 

IV.  Zero  Scan  Length 

In  the  idealized  case  of  the  zero  scan  length,  all  reports  from  a  given  scan  have  identical 
timestamps.  To  perform  gating,  the  track  file  is  projected  to  the  time  of  the  reports, 
then  the  search  radius  will  not  depend  on  object  dynamics  since  =  0  in  Eq.  (3).  For 
the  solution  of  the  equal  timestamp  case,  therefore,  we  calculate  a  seeo'ch  radius  from  a 
given  score  threshold.  More  specifically,  in  doing  a  search  per  scan  per  report  j  on  a  track 
file,  we  will  determine  the  following:  given  a  score  threshold  Sj,  calculate  the  radius 
per  report,  Ro  =  Ro  =  Ro{Sj),  such  that  all  the  report-track  pairs  which  are  separated 
by  distances  larger  than  Rq  will  not  have  scores  above  the  threshold.  Calculation  of 
this  radius  is  complicated  by  the  fact  that  the  tracks  in  our  database  have  nonidentical 
covariance  matrices,  or  hyperellipses,  of  nontrivial  distribution  in  volume  and  shape. 
Consider,  however,  the  case  where  d  =  1.  For  a  one  dimensional  Gaussian  function 
with  a  fixed  F,  the  function  value  is  strictly  nonincreasing  as  its  argument  R?  increases 
and  falls  below  the  threshold  after  R?  gets  laiger  than  some  value,  say  Rq.  One  can 
then  determine  analytically  an  Rq  that  is  independent  of  the  F  distribution.  Figure  1 
below  is  a  contour  plot  for  the  function  5(7, r)  =  exp(— r^/27)/\/(27r7).  The  lines  in 
the  figure  are  actually  iso-score  lines.  Notice  that  each  exhibits  a  maximum  radius  r,nax 
such  that  for  scores  exceeding  the  threshold,  all  possible  r  are  smaller  than  Tmax-  This 
rmoi  is  a  suitable  choice  for  Rq.  Appendix  I  has  a  derivation  for  a  useful  search  radius 
for  dimension  d  >  1;  the  result  is; 
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where  (aj)  is  an  upper  bound  on  the  maximum  eccentricity  among  the  error  ellipses  in 
a  scan. 

Equation  (A12)  gives  a  suitable  choice  for  Rq  since;  (1)  all  tracks  which  can  correlate 
with  observation  j  will  be  within  Rq  of  the  observation,  and  (2)  it  is  a  track-independent 
radius  and  thus  can  be  used  for  a  search  on  a  track  file.  We  should  note  that  other  pre¬ 
scriptions  for  calculating  time  independent  gate  volumes  might  meet  the  above  criteria. 
Such  prescriptions  might  differ  in  how  close  Vc{Ro)  approaches  Vgi  and  how  expensive 
it  will  be  to  calculate  Vg{Rq).  Preliminary  analysis  indicates  that  our  Vg{Rq)  is  a  rather 
coarse  gate  that  is  inexpensive  to  calculate. 

Though  the  equal  timestamp  behavior  of  the  search  algorithm  is  straightforward  enough, 
we  include  Plot  1  with  data  from  a  simulation  (see  Appendix  II).  Plot  1  verifies  the  search 
time  scaling  for  an  unequal  timestamp  case  which  is  not  quite  the  r  =  0  limit,  but  with 
relatively  small  differences  in  timestamp.  In  this  case  r  and  our  overestimate  of  Via 
were  small  enough  so  that  the  dominant  term  in  the  scaling  of  the  BLD  tree  is  seen 
to  be  close  to  the  0{NR\ogNr)  diagonal  line  in  the  figure.  The  (9(A^r log A^r)  CPU 
expense  of  creating  the  data  structure  was  also  verified. 


V.  The  General  Algorithm  and  Analysis  for  the  u|^r|  >>  Rq  Case 

When  the  reports  have  unequal  timestamps  what  are  the  best  times  to  which  to  integrate 
the  tracks  for  making  the  potential  pair  matching?  Following  the  textbook  procedure 
of  making  all  the  tracks  valid  at  the  beginning  or  all  valid  at  the  end  of  the  scan  can 
result  in  the  case  of  Eq.  (4b)  and  possibly  in  a  combinatorial  bottleneck.  We  have  made 
the  following  observations:  (1)  some  search  structures  are  cheap  to  make,  CPU  wise, 
and  (2)  by  making  various  copies  of  the  track  data  structure  -  the  various  copies  valid 
at  different  times  -  a  tradeoff  can  be  made  between  the  time  spent  on  the  creation  of 
the  data  structures  and  the  time  spent  on  searching  and  scoring.  We  shall  see  that  the 
payoff  from  spending  more  time  in  creating  data  structures  increases  as  the  scan  length 
contribution  to  Vg  increases,  given  that  the  other  relevant  parameters  are  held  fixed. 

If  we  integrate  in  time  to  make  Md  sequential  track  data  structures  (TDS)  valid  at 
Md  equally  spaced  times  within  the  scan  of  length  t  (See  Fig.  2),  then  any  report 
would  be  at  most  |^T(  =  tI{2Md)  time  units  away  from  a  TDS.  The  average  radius 
for  the  search  is  then  decreased  by  the  factor  Mu  as  compared  to  the  case  where  we 
have  one  TDS  copy  at  the  middle  of  the  scan.  Therefore  the  volume  extent  as  well  as 
the  average  number  of  candidates  returned  {Nq)  is  smaller  by  (1/Md)*  in  the  isotropic 
dense  limit  /-dimensional  case.  More  precisely,  assume  that  the  density  p  of  objects 
in  space  is  constant  and  uniform.  Then  the  average  number  of  candidate  tracks  found 
for  each  report  depends  on  the  average  search  volume  Vg  =  7(^)(^g)’  where  7(f)  is  a 
geometric  factor  depending  on  the  dimension  (  of  the  report  state  vector.  The  brackets 
(•)  denote  the  average  over  the  temporal  range  (/,  —  ST,  /,  -f  ST),  where  i{R)  labels  the 
TDS  selected  to  be  closest  to  the  report.  Assuming  that  the  time  distribution  of  reports 
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Since  the  searching  time  and  the  scoring  time  depend  on  Nc,  the  scoring  time  being 
directly  proportional  to  No,  the  use  of  multiple  extrapolated  track  files  (Md  >  1)  to 
cover  the  scan  interval  can  reduce  the  cost  of  the  gating  process  by  reducing  the  search 
volume. 

Plot  2  shows  an  example  of  poor  scaling  due  to  a  “large”  scan  length.  The  trajectories 
which  generated  the  data  for  Plot  2  are  equivalent  to  those  for  Plot  1  except  that  the 
scan  length  is  increased  by  a  factor  of  five  and  therefore  the  average  search  radius  is  also 
increased.  In  this  case  the  search  algorithm’s  visible  CPU  cost  was  not  OiNulogNr)- 
For  Nt  =  Nn  =  32A"',  the  search  time  increased  by  a  factor  5,  and  the  number  of  near 
neighbors  returned  was  increzised  by  a  factor  of  about  10.  Plot  3  shows  how  the  scaling 
of  the  search  time  behaves  for  the  same  report  scenario  as  in  Plot  2,  except  that  now 
Md  =  5  for  the  datasets  of  4K  and  32K  tracks.  The  best  fit  lines  approach  the  log  Nt 
scaling  line  more  closely.  For  the  same  runs  and  Nj  =  Njt  =  32/v ,  the  number  of  near 
neighbors  returned,  and  therefore  the  number  of  pairs  which  would  be  evaluated  with 
Eq.  (1),  decreased  by  an  average  factor  of  9.8  when  Md  was  changed  from  1  to  5.  This 
reduction  also  signifies  that  the  number  of  tracks  per  report  found  to  meet  the  time 
dependent  gate  volume  was  nearly  reduced  to  its  simulation  optimum  average  value  of 
1.  This  factor  of  9.8  was  found  to  be  62.7  if  the  one  TDS  was  placed  at  the  time  of  the 
beginning  of  the  scan  instead  of  the  recommended  placement  for  Md  —  l! 

Was  the  price  paid  for  overhead  of  the  multiple  data  copies  too  high  so  that  the  overall 
cost  would  still  scale  poorly?  Let  Nq  be  the  report  averaged  number  of  near  neighbors 
returned.  Then  the  total  CPU  cost  can  be  modeled  as 


C{Nt,N  r.MdiN  a)  —  CeMDNT  +  CdMDNT\ogNT  +  CteNii{\ogNT  +  N  o)  +  CgcNGNfi 

(9) 

where  the  terms  on  the  RHS  of  Eq.  (9)  give,  respectively,  the  cost  for  integrating  the 
tracks  to  the  desired  time  of  the  data  structures,  the  cost  of  making  the  tree  data 

^An  overestimate  of  the  worst  case  is  when  the  reports  are  8T  =  t/{2Md)  time  units  away  from  a 
TDS,  i.e.  at  the  furthest  possible  time  difference,  where  the  reported  average  case  is  small  by  2^/(^+  I) 
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structures,  the  cost  of  searching  the  appropriate  tree  data  structure  for  each  report,  and 
the  cost  of  scoring  the  pairs. 

Equation  (9)  with  our  model  for  Eq.  (7)  has  a  minimum  value  which  occurs  for  the 
following  optimal  Md 


Mdo  =  (10) 

Kl 

where 

Ki  =  NT{C,  +  Cj\ogNT), 

=  +  .  (11) 

and  the  total  cost  is 

Ctnin  ~  CaeA^filog  AV  +  Mgo  ~^~(C'eA^r  +  CrfA^r  log  A’t)-  (12) 

The  above  equation  depends  on  the  important  parameters  for  scaling  of  multitar^t 
tracking,  as  previously  argued,  except  that  it  does  not  depend  on  combinations  of  Rq 
with  or  because  of  the  approximation  in  Eq.  (4b). 

VI.  Analysis  of  the  General  Algorithm 

Let  ^  =  ar I{2RqMd),  where  the  symbol  for  the  number  of  TDS  is  now  to  make 
a  distinction  for  the  limit  of  Eq.  (4b).  Instead  of  taking  the  approximation  in  Eq.  (6a) 
leading  to  Eq.  (6b),  use  the  binomial  expansion  and  Nq  =  pVa  to  obtain  from  Eq.  (6a): 

^  ^  i=i 

To  find  Cmin  it  is  useful  to  find  the  partial  of  Eq.  (13)  with  respect  to  Ad/?: 


dNcjdAio  = 


7(^)  {RqY 

+  1  Md 


(14) 


This  allows  us  to  obtain  an  equation  for  the  number  of  TDS  that  minimize  the  cost  of 
Eq.  (9): 


(15) 


The  above  is  a  polynomial  equation  for  Ad do  of  £  terms  and  of  degree  £  +  1.  For  1—1, 

Ad DO  =  Mdo- 

We  now  verify  the  Cmin  of  Eq.  (1^)  (evaluated  when  or  >>  2^)  is  an  overestimate 
when  it  is  not  true  that  or  >>  TRp.  We  compare  Adpp  and  Mpo  for  a  given  Nj  and 
Nr,  and  with  a  fixed  estimate  of  pV ic  through  pVq{Ro).  Let  be  some  value  of  a 
scan  length  for  which  ar^  >>  2Ro  and  for  which  Eq.  (10)  was  evaluated  to  be  Mooiri). 
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Then,  using  the  ratio  of  Eq.  (15)  at  r  and  at  tl  and  using  Eq.  (10)  for  Aipo  at  r/,,  the 
value  of  Aino  at  some  arbitrary  scan  length  is: 

Algo'  E de) 

TL  t  or 

where  the  magnitude  in  the  approximation,  due  to  MdoI'^'l)  =  ^^Doi'^L)  is  given  by 
Eq.  (15).  The  equation  above  shows  how  to  calculate  Mdo  for  an  arbitrary  scan  length 
given  that  it  is  evaluated  for  r^.  Also,  Eqs.  (9),  (13)  and  (16)  give  the  cost  of  the 
gating  process  in  terms  of  relevant  parameters  Nx,  p,  r,  Ro,  (and  therefore  Ro) 
and  combinations  of  them.  Notice  that  each  of  the  terms  in  Eqs.  (13)  and  (16)  have 
their  contribution  in  r  in  the  form  of  r*,  where  i  is  some  positive  integer.  Thus  No  (and 
Ng)  and  decrease  as  r  decreases  on  some  interval  (0,tx,).  Notice  also  that  as  No 

and  Mdo  decrease,  the  cost  as  given  by  Eq.  (9)  decreases.  And  since  for  t  =  ti,  Mdo 
=  Mi5o(l  +  higher  orders  in  2Ro/a'rL),  the  large  scan  length  case  cost  is  an  overestimate 
of  the  cost  for  the  general  case  with  a  smaller  scan  length  and  with  the  other  parameters 
held  fixed. 

VII.  Discussion 

Our  choice  of  the  time  and  dynamics  dependent  component  of  Ro  in  Eq.  (3)  might 
not  be  the  only  choice  possible  in  that  other  models  might  exist  which  allow  a  range 
search.  However,  making  multiple  projections  of  the  data  to  be  searched  seems  to  be  an 
option  for  consideration  whenever  the  bounds  on  the  search  range  are  large  due  to  large 
timestamp  differences.  Related  to  this  is  that  not  all  of  the  search  algorithms  require 
that  the  coordinate  ranges  per  search  be  identical  in  all  directions.  It  might  be  possible 
to  separate,  by  directional  components,  both  Rq  and  the  time  dependent  component  of 
Rq.  What  is  required  is  that  one  search  range  per  coordinate  -  and  per  report  when 
doing  a  search  on  tracks  -  be  given  to  the  search  algorithm  such  that  all  pairs  which 
meet  Eq.  (1)  are  satisfied.  If  it  is  not  true  that  on  the  average  pVciRo)  ^  where  N 
is  the  number  of  tracks  when  doing  a  search  on  tracks  and  where  Vg{Ro)  is  determined 
as  in  our  paper,  it  might  be  profitable  to  separate  Vg(Ro)  into  search  components  or 
to  calculate  Vg{Ro]  by  a  different  prescription.  However,  it  should  be  noted  that  Eqs 
(10)  and  (12)  are  independent  of  Rq  and  so  the  previous  suggestion  might  not  improve 
execution  speed  in  the  large  scanlength  case.  On  the  other  hand,  execution  speed  in  the 
large  scanlangth  case  might  be  noticably  improved  by  a  search  space  model  that  will 
allow  one  to  seperate  Vg{ST).  Note  that  an  algorithm  that  uses  velocity  information 
may  reduce  the  search  volume  to  a  half-sphere  or  p-irhaps  even  a  cone  or  a  rectangle. 
This  has  the  effects  of  changing  7(f)  but  leaving  Eq.  (12),  and  therefore  scaling,  as  is. 

In  addition  to  the  above  considerations,  the  value  of  the  optimal  number  for  Md  depends 
on  the  computer  used,  and  relative  software  component  speeds,  through  the  constants 
in  Eq.  (9).  For  example  a  vector  computer  with  good  floating  point  hardware  will 
no  doubt  reduce  C41  over  a  generic  computer;  but  because  Ci  can  be  reduced  by  an 
even  larger  factor  since  vectorization  of  the  integration  of  the  equations  of  motion  will 
be  straightforward  and  possibly  complete.  This  will  incline  one  to  have  a  larger  Mp. 
Performance  discussion  of  different  hardware/software  setup  can  get  quite  detailed  and 
possibly  in-house  experimentation  will  be  necessary  for  determining  an  optimal  Mp. 
Finally,  it  should  be  noted  that  Eq.  (9)  was  not  presented  with  the  obvious  constraint 
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that  Md  he  an  integer  greater  than  zero.  Rounding  up  to  an  integer  is  advisable  in 
implementation  when  using  Eqs.  (10)  or  (16). 

VIII.  Conclusion 

In  this  paper  we  have  presented  an  efficient  approach  to  gating  in  multiple-target  track¬ 
ing.  It  is  general  in  the  sense  that  it  will  perform  with  good  scaling  even  for  the  case  of 
unequal  observation  timestamps.  An  analysis  has  been  provided  detailing  how  the  al¬ 
gorithm  scales  as  a  function  of  the  important  parameters:  AV,  A'h,  r,p  and  an  estimate 
of  the  time  independent  gate,  Vc{Rq).  The  algorithms  has  a  worst  case  scaling,  in  the 
timestamp  parameter,  as  follows: 

C'min  =  0{ Nr  -f  MoqNt  log  A'r ), 

where  Mdq  «  {pN rt^  j Nt{c ■{■log  At))'+>  .  It  is  also  shown  how  to  calculate  an  optimum 
number  of  search  data  structure  in  general,  and  how  the  cilgorithm  scales  in  generail. 

Appendix  I  :  Derivation  of  Search  Radius 
Covariances  and  Score  Threshold  Contribution 

Given  a  report  at  observation  time  t  and  a  set  of  tracks  which  have  been  projected  to 
time  t,  the  gating  process  provides  the  subset  of  tracks  whose  scores  exceed  a  threshold 
value,  denoted  here  as  5nun.  The  successful  application  of  efficient  search  algorithms 
to  improve  scaling  with  the  number  of  objects  requires  that  a  search  region  (or  search 
radius  Rq)  be  defined  independent  of  the  particular  track  when  doing  a  search  upon 
the  track  database.  A  satisfactory  expression  for  the  search  radius,  for  example,  must 
ensure  that  all  acceptable  tracks  fall  within  a  distance  Rq  of  the  report.  This  appendix 
derives  such  an  equation  for  the  search  region. 

Consider  a  given  report  and  track  pair  whose  respective  labels  are  i  and  j,  whose  M- 
dimensional  measurement-prediction  difference  is  denoted  in  Eq.  (1),  and 

whose  position  displacement  is  denoted  by  6R  =  5R,j  .  To  obtain  a  search  criterion, 
use  the  axiom  that  if  the  score  Sij  >  5nuni  then  the  maximum  possible  score  So  (given 
^X)  must  also  exceed  5inin-  For  this  value  of  6X.,  therefore,  we  compute  the  covariance 
matrix  Fq  which  maximizes  Si^  The  matrix  Fq  depends  directly  on  ^X,  and,  therefore, 
so  does  Sq.  Comparison  of  Sq  with  5min  establishes  a  necessary  condition  between  6X 
and  Smin  for  all  satisfactory  track-report  pairs.  This  condition  gives  an  upper  bound 
on  the  search  area  or  radius,  but  not  necessarily  a  least  upper  bound.  For  each  report, 
one  can  identify  a  subset  of  tracks  which  should  contain  all  tracks  meeting  the  scoring 
criteria  as  well  as  a  small  number  which  do  not.  Performing  this  identification  for  all 
reports  results  in  a  set  of  candidate  track-report  pairs.  Computing  the  actual  scores  of 
these  pairs  and  selecting  those  whose  scores  exceed  the  threshold  5min  completes  the 
gating  process. 

We  begin  with  Eq.  (1)  and  notice  that  both  the  argument  of  the  exponential  in  the 
numerator  and  the  determinant  of  the  the  covariance  matrix  F  in  the  denominator  are 
invariant  under  similarity  transformations  of  the  coordinates.  Assume  that  a  transfor¬ 
mation  exists  which  renders  F  diagonal  with  nonzero  elements  given  by  the  eigenvalues 
Aj,  A2, . . . ,  A\f .  Here  d  is  the  dimension  of  the  (Euclidean)  report  state  vector  space. 
Then  for  a  given  report-track  pair,  the  score  is 
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-log5  =  r^-7^  +  -log2jr  +  -^logA,  .  (yll) 

Note  that  the  vector  (JXhas  also  undergone  transformation,  although  we  retain  the 

same  notation  for  convenience.  For  a  given  value  of  <5X,  we  now  wish  to  find  the  set 
of  eigenvcdues  {A*}  which  maximize  the  score.  This  requires  that  w’e  take  the  partial 
derivative  of  the  right  hand  side  of  Eq.  (Al)  with  respect  to  each  eigenvzJue  A*  and  set 
the  resulting  expression  equal  to  zero.  In  this  way  we  obtain  the  condition 

Xk  =  SXl  Wk  .  {A2) 


Substituting  this  in  Eq.  (Al)  gives  us 

-logSo=^(l  +  log27r+-X:iog6At)  .  (A3) 

From  Eq.  (A3)  and  the  necessary  condition  for  a  given  track-report  pair  to  have  a  score 
exceeding  5mini  we  obtain  the  following  condition  for  the  search  volume  (region): 

(«?>= (/IS) 

fc=i 


is  the  geometric  average  of  the  components  of  the  difference  between  the  track  and 
report  positions  in  the  basis  which  diagonalizes  the  covariance  matrix. 

Equations  (A4-5)  constitute  our  most  general  result.  This  caise,  is  not  straightforward 
in  allowing  the  definition  of  a  search  radius,  since  the  geometric  average  and  arithmetic 
average  of  the  squared  components  of  the  vector  ^X  satisfy 

^  (.6) 

by  a  theorem  in  mathematical  analysis.  One  can  see  the  reason  by  considering  the  cir¬ 
cumstance  in  which  the  sensors  have  isotropic  resolution.  Then  the  covariance  matrices 
will  have  d-fold  degenerate  eigenvalues  so  that 

Afc  =  -^I^Xj^  (isotropic  case)  (A7) 

a 

will  maximize  the  score  in  Eq.  (Al)  and  the  result  corresponding  to  Eq.  (A4)  is 
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I^RI"  <  1^X1^  <  ^(5™n)--  =  Rl  (^8) 

Because  the  magnitude  of  the  track-report  displacement  appears  on  the  left-hand  side 
(LHS),  the  expression  on  the  right-hand  side  (RHS)  constitutes  a  search  radius,  denoted 
as  Rq. 


For  the  more  general  case,  it  Is  necessary  to  have  an  estimate  of  the  limits  of  anisotropy 
of  the  set  of  residual  covariances  for  the  coordinate  system  where  the  near  neighbors 
search  is  done.  Overestimates  can  be  obtained  if  the  set  of  eigenvalues  for  all  the  residual 
covariances  are  known.  Obtaining  the  eigenvalues  is  an  0(Nt)  operation  which  for  two 
and  three  dimensions  does  not  have  to  be  expensive,  as  we  will  suggest  latter,  (this 
shoiild  be  a  note.)  For  the  jth  track  in  the  databzise,  let  the  ratio  of  the  fcth  principal 
axis  to  the  longest  principal  axis  be  given  by  where  Ok^  <  1.  The  eigenvalues  of 
matrix  Fj  should  then  satisfy 


Xi  =  af  V  (A9) 

where  A-’  is  the  largest  eigenvalue  of  Fj.  We  now  choose  the  Fj  which  has  the  smallest 
geometric  mean  of  the  s.  For  this  F,  then,  with  Eq.  (A9)  and  Eq.  (Al)  and  computing 
the  optimum  value  of  A  gives  us 


=  •  (^n) 

Denoting  the  geometric  average  as  before,  noting  the  necessary  condition  on  5o,  and 
remembering  that  <  1,  we  obtain 

Again  since  |^R|*  appears  on  the  LHS,  the  extreme  RHS  defines  a  search  radius  Rq. 

Appendix  II  :  Scenario  Description 

For  the  sake  of  completeness  we  state  that  the  basic  scenario  presented  in  this  paper 
consisted  of  Nt  =  Np,  reports  and  tracks  obtained  by  adding  Gaussian  noise  to  ground 
truth  trajectories.  The  ground  truth  was  in  an  XYZ  common  coordinate  system.  The 
trajectories  were  Newtonian  and  paraboDc  with  a  constant  gravitational  force  in  the  Z- 
direction;  also  they  were  determined  by  choosing  initial  positions  and  final  destinations 
random  in  X,Y  with  Z  =  0  and  bounded  on  the  X  —  Y  plane  with  one  rectangle  for 
the  initial  positions  and  another  for  the  final  destinations.  The  two  rectangles  were  set 
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sufficiently  far  apart  so  that  the  average  flight  would  be  about  25  minutes  in  duration 
with  the  ground  truth  objects  flying  at  about  11  km/sec.  A  scan  consisted  of  a  sampling 
on  the  trajectories  in  the  form  of  reports  with  timestamps  within  a  scan  of  length  r. 
The  only  parameters  changed  to  generate  all  the  data  presented  in  the  plots  were  Np, 
(and  Nt),  t  iMd  (the  number  of  track  data  structure  copies),  and  the  object  density 
as  a  consequence  of  changing  Np  (and  Nj)- 
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Fig.  1  —  Contour  plot  of  the  score  function  s  =  e  2iro .  The  Tq's  jre  the  chosen  suitable  time  independent  radius 

component  due  to  covariances  and  score  threshold  contributions. 
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